SUBROUTINE G_ENE
!-----------------------------------------------------------------------
!
! Definition von G:  Energiegleichung
!

      use primvar,  only : G, X, XZ, XA, MR, MD, ME, MUr, MJ, XZ, Gnorm
      use physco,   only : sqrt2pi, pi, z12, z2, z4, z8, StBolz
      use global,   only : tst, relax
      use config,   only : np, rkind
      use zvar,     only : H_pZ, H_pZ_flag
      use geomvar,  only : S_flux, S_volZ, S_vol, S_volA
      use matvar,   only : Pgas0Z, OPAplkZ, TgasZ
      use viscvar,  only : muQ_r, Uterm_r, muQ_p_ene, Uterm_p
      use advecvar, only : e_adv
      use eddvar,   only : mu_z
      use RBvar,    only : L_star, H_star, f_irr, e_in, e_out


      implicit none

      integer       :: i
      logical, save :: init = .true.

      if (init) then
         if (.not. relax) then
            write(66,"(a)") "G_ene.f90:   Innenrand: de/dR = 0 ;  Aussenrand: de/dR = 0"
         else
            write(66,"(a)") "G_ene.f90:   Zum Relaxieren: Innenrand: e = e_in ;  Aussenrand: e = e_out"
         end if
         init = .false.
      end if


!-----------------------------------------------------------------------
!    Randbedingungen
!-----------------------------------------------------------------------

      ! innere Pseudozellen: i=1 & i=2
!~       if ( X(MUr,3) * tst / X(MR,3) > +1.e-30_rkind ) then  ! Einstroemung ueber Innenrand

         G(ME,1)        =  X(ME,1)    - e_in
         Gnorm(ME,1)    =  X(ME,1)
         G(ME,2)        =  X(ME,2)    - e_in
         Gnorm(ME,2)    =  X(ME,2)

!~       else
!~
!~          G(ME,1)        =  X(ME,1)    - X(ME,2)
!~          Gnorm(ME,1)    =  X(ME,1)
!~          G(ME,2)        =  X(ME,2)    - X(ME,3)
!~          Gnorm(ME,2)    =  X(ME,2)
!~
!~          e_in           =  X(ME,3)
!~
!~       end if


      ! aeussere Pseudozellen + ueberzaehlige skalare Zelle: i=np, i=np-1, i=np-2
!~       if ( X(MUr,np-2) * tst / X(MR,np-2) < -1.e-30_rkind ) then  ! Einstroemung ueber Aussenrand

         G(ME,np)       =  X(ME,np)   - e_out
         Gnorm(ME,np)   =  X(ME,np)
         G(ME,np-1)     =  X(ME,np-1) - e_out
         Gnorm(ME,np-1) =  X(ME,np-1)
         G(ME,np-2)     =  X(ME,np-2) - e_out
         Gnorm(ME,np-2) =  X(ME,np-2)

         ! Wenn man e_out als RB setzt wird der Einstrahlungsterm in der aeussersten skalaren Zelle (idx = np-3) potentiell unphysikalisch.
         ! Mit folgender RB fuer idx = np-3 kann man dieses Problem umgehen.
         G(ME,np-3)     =  X(ME,np-3) - X(ME,np-4)
         Gnorm(ME,np-3) =  X(ME,np-3)

!~       else
!~
!~          G(ME,np)       =  X(ME,np)   - X(ME,np-1)
!~          Gnorm(ME,np)   =  X(ME,np)
!~          G(ME,np-1)     =  X(ME,np-1) - X(ME,np-2)
!~          Gnorm(ME,np-1) =  X(ME,np-1)
!~          G(ME,np-2)     =  X(ME,np-2) - X(ME,np-3)
!~          Gnorm(ME,np-2) =  X(ME,np-2)
!~
!~          e_out          =  X(ME,np-3)
!~
!~       end if


!-----------------------------------------------------------------------
!    Restlicher Bereich
!-----------------------------------------------------------------------

      do i=3,np-4

         G(ME,i) = S_vol(i) * X(ME,i) * X(MD,i) - S_volA(i) * XA(ME,i) * XA(MD,i)                                    & ! temporal diff.
                 + S_flux(i+1) * e_adv(i+1) - S_flux(i) * e_adv(i)                                                   & ! advection
                 + sqrt2pi * H_pZ(i) * Pgas0Z(i) * z2*pi * ( XZ(MR,i+1)*XZ(MUr,i+1) - XZ(MR,i)*XZ(MUr,i) ) * tst     & ! Pgas div(u) term
                 - z4*pi * OPAplkZ(i) * XZ(MD,i) * ( XZ(MJ,i) - StBolz/pi * TgasZ(i)**4 ) * S_volZ(i) * tst          & ! (J-S) term
                 - z12 * ( muQ_r(i) * Uterm_r(i)**2 + muQ_p_ene(i) * Uterm_p(i)**2 ) * S_volZ(i) * tst               & ! viscosity
                 - L_star * f_irr * H_pZ_flag(i) * (   (H_pZ(i+1) - H_star) / XZ(MR,i+1)                             & ! irradiation
                                                     - (H_pZ(i)   - H_star) / XZ(MR,i)   ) * tst                     & ! ...
                 + z8*pi * mu_z * XZ(MJ,i) * S_volZ(i) * tst                                                           ! radiative cooling


         Gnorm(ME,i) = max(                                                                                          &
              abs( S_vol(i) * X(ME,i) * X(MD,i) ),                                                                   & ! cell content
              abs( S_flux(i+1) * e_adv(i+1) ), abs( S_flux(i) * e_adv(i) ),                                          & ! advection
              abs( sqrt2pi * H_pZ(i) * Pgas0Z(i) * z2*pi * ( XZ(MR,i+1)*XZ(MUr,i+1) - XZ(MR,i)*XZ(MUr,i) ) * tst ),  & ! Pgas div(u) term
              abs( z4*pi * OPAplkZ(i) * XZ(MD,i) * max( XZ(MJ,i) , StBolz/pi * TgasZ(i)**4 ) * S_volZ(i) * tst ),    & ! (J-S) term
              abs( z12 * ( muQ_r(i) * Uterm_r(i)**2 + muQ_p_ene(i) * Uterm_p(i)**2 ) * S_volZ(i) * tst ),            & ! viscosity
              abs( L_star * f_irr * H_pZ_flag(i) * (   (H_pZ(i+1) - H_star) / XZ(MR,i+1)                             & ! irradiation
                                                     - (H_pZ(i)   - H_star) / XZ(MR,i)   ) * tst ),                  & ! ...
              abs( z8*pi * mu_z * XZ(MJ,i) * S_volZ(i) * tst )                                                       & ! radiative cooling
                          )

! In der Strahlungsenergiegleichung stellt J-S die wesentlichste Bedingung zur Bestimmung von J dar,
! im groessten Teil des Sterns wird J direkt durch diese Bedingung definiert - alle dynamischen terme sind dort vernachlaessigbar.
! Deshalb ist es notwendig den J-S Term fuer die Normierung der Strahlungsenergiegleichung aufzuspalten.
!
! Im Gegensatz dazu wirkt J-S hier in der Energiegleichung nur als Quellterm - man koennte ihn also als Ganzes fuer
! die Normierung nehmen.
! Allerdings kann das fuer stationaere Loesungen, d.h. u~0, zu Problemen fueren. In diesem Fall fallen fast alle anderen
! Terme weg und uebrig bleibt nur noch    e - eA = -4 pi kappa (J-S) tst
! Da J-S nur mit begrenzter Genauigkeit bekannt ist wird diese Gleichung fuer grosse Zeitschritte unloesbar, die man aber
! natuerlich gerade fuer stationaere Loesungen anstrebt, z.b. beim Relaxieren.
! Im Prinzip ist das genau das gleiche Problem das uns G_ene_radene gefuehrt hat - nur das es hier schon fruher auftritt.
! Folglich habe ich auch hier den J-S term fuer die Normierung aufgespalten.

      end do


END SUBROUTINE G_ENE

